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Abstract 

This  two-part  paper  discusses  robustification  methodologies  for  linear-iterative  distributed  algorithms  for 
consensus  and  coordination  problems  in  multicomponent  systems,  in  which  unreliable  communication  links  may 
drop  packets.  We  consider  a  setup  where  communication  links  between  components  can  be  asymmetric  (i.e., 
component  j  might  be  able  to  send  information  to  component  i,  but  not  necessarily  vice-versa),  so  that  the 
information  exchange  between  components  in  the  system  is  in  general  described  by  a  directed  graph  that  is 
assumed  to  be  strongly  connected.  In  the  absence  of  communication  link  failures,  each  component  i  maintains 
two  auxiliary  variables  and  updates  each  of  their  values  to  be  a  linear  combination  of  their  corresponding 
previous  values  and  the  corresponding  previous  values  of  neighboring  components  (i.e.,  components  that  send 
information  to  node  i).  By  appropriately  initializing  these  two  (decoupled)  iterations,  the  system  components 
can  asymptotically  calculate  variables  of  interest  in  a  distributed  fashion;  in  particular,  the  average  of  the  initial 
conditions  can  be  calculated  as  a  function  that  involves  the  ratio  of  these  two  auxiliary  variables.  The  focus 
of  this  paper  to  robustify  this  double-iteration  algorithm  against  communication  link  failures.  We  achieve  this 
by  modifying  the  double-iteration  algorithm  (by  introducing  some  additional  auxiliary  variables)  and  prove  that 
the  modified  double-iteration  converges  almost  surely  to  average  consensus.  In  the  first  part  of  the  paper,  we 
study  the  first  and  second  moments  of  the  two  iterations,  and  use  them  to  establish  convergence,  and  illustrate 
the  performance  of  the  algorithm  with  several  numerical  examples.  In  the  second  part,  in  order  to  establish  the 
convergence  of  the  algorithm,  we  use  coefficients  of  ergodicity  commonly  used  in  analyzing  inhomogeneous 
Markov  chains. 
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I.  Introduction 

The  design  of  protocols  and  algorithms  for  distributed  computation  and  control/decision  tasks  has 
attracted  significant  attention  by  the  computer  science,  communication,  and  control  communities  (e.g., 
[1],  [2],  [3],  [4],  [5],  [6],  [7]  and  references  therein).  For  example,  given  i)  a  collection  of  robots 
moving  in  the  plane,  ii)  a  collection  of  sensors  in  a  sensor  network,  or  iii)  a  collection  of  distributed 
energy  resources  in  an  electrical  grid,  the  components  may  be  interested  in,  respectively,  i)  agreeing 
on  a  common  direction  to  follow  (this  common  direction  could  be  provided  by  a  leader  robot),  ii) 
measurement  averaging  (with  each  sensor  providing  a  local  measurement  of  a  global  quantity),  or  iii) 
collectively  providing  a  predetermined  total  amount  of  active  power  subject  to  the  constraints  of  each 
distributed  resource.  In  the  control  literature,  the  first  and  second  problems  are  respectively  known  as 
consensus  and  average  consensus  (see,  e.g.,  [2]),  whereas  the  third  problem  can  be  considered  as  a 
distributed  resource  coordination  problem  [8],  [9]. 

In  this  two-part  paper,  we  consider  multicomponent  systems  in  which  each  component  can  exchange 
information  with  other  components  in  its  neighborhood  in  order  to  compute,  in  a  distributed  fashion, 
some  quantity  of  interest.  In  our  setup,  communication  links  between  components  (nodes)  can  be 
asymmetric  (i.e.,  component  j  might  be  able  to  send  information  to  component  i,  but  not  necessarily 
vice-versa),  a  situation  that  arises  in  a  wireless  setting  if  the  transmission  power  available  to  different 
nodes  are  also  different.  In  this  setting,  the  information  exchange  between  components  in  the  system 
can  be  described  by  a  directed  graph  which  is  assumed  to  be  strongly  connected.  Through  an  iterative 
process,  nodes  in  the  network  are  required  to  compute  (using  only  information  made  available  by  their 
neighbors)  the  quantity  of  interest.  In  particular,  we  study  linear-iterative  algorithms  in  which  each  node 
j  maintains  a  value  (or  a  set  of  values)  that  is  updated  to  be  a  weighted  linear  combination  of  node  j’s 
own  previous  value  and  the  previous  values  of  its  neighboring  nodes  (i.e.,  nodes  that  transmit  information 
to  node  j).  The  main  focus  of  the  paper  is  to  develop  strategies  to  robustify  the  linear-iterative  algorithms 
described  above  against  communication  links  that  may  drop  packets. 

In  the  context  of  consensus  and  average-consensus  problems,  an  extensive  literature  in  the  control 
community  focuses  on  the  linear-iterative  algorithms  described  above  (e.g.,  [7],  [10],  [2],  [11],  [12], 
[13],  [14]  and  references  therein).  These  works  have  revealed  that  if  the  network  topology  satisfies 
certain  conditions,  the  weights  for  the  linear  iteration  can  be  chosen  so  that  all  the  nodes  asymptotically 
converge  to  the  same  value  (even  if  the  network  connections  are  time- varying).  Additionally,  if  the 
interconnection  topology  is  invariant  and  bidirectional  (i.e.,  if  node  j  can  send  information  to  node  i, 
then  node  i  can  send  information  to  node  j),  simple  techniques  can  be  used  to  choose  the  weights  of  the 
linear  iteration  so  as  to  ensure  that,  after  running  the  linear  iteration,  the  nodes  will  asymptotically  reach 
consensus  to  the  average  of  their  initial  values  [2],  [11],  [12].  Other  works  have  looked  at  the  consensus 
and  average-consensus  problems  when  the  interconnection  topology  is  described  by  a  directed  graph. 
In  particular,  the  authors  of  [15]  focus  on  continuous-time  linear  iterations  and  state  necessary  and 
sufficient  conditions  for  a  network  of  integrators  to  asymptotically  reach  agreement  to  a  common  value 
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(but  not  necessarily  the  average  of  their  initial  values).  Similarly,  the  authors  of  [12]  consider  discrete¬ 
time  iterations,  and  provide  necessary  and  sufficient  conditions  on  the  weights  that  allow  the  nodes  to 
asymptotically  reach  consensus  to  the  average  of  their  initial  values.  Additionally,  the  work  in  [16], 
[17]  discusses  how  average-consensus  can  be  reached  asymptotically  with  linear- iterative  algorithms  in 
which  the  nodes  use  fixed  weights  in  their  linear  updates  and  also  develops  linear-iterative  algorithms 
where  the  nodes  adapt  their  weights  in  a  distributed  fashion  so  that  asymptotically  average-consensus 
is  reached.  In  the  context  of  resource  coordination,  there  is  some  recent  work  [18],  [8],  [9]  that  also 
focuses  on  linear-iterative  algorithms,  similar  to  those  used  to  address  consensus  and  average-consensus 
problems.  Other  recent  work  has  addressed  the  related  problem  of  achieving  consensus  and  average- 
consensus  in  a  multicomponent  system  where  some  nodes  can  exhibit  malicious  behavior  [19],  [20]. 
These  works  assume  fault-free  communication  links,  but  are  related  to  what  we  do  in  this  paper  in  the 
sense  that  they  can  be  used  to  handle  unreliable  nodes  (as  opposed  to  links). 

In  our  development,  we  adopt  a  very  general  model  for  the  communication  modality  between  nodes, 
which  allows  asymmetric  information  structures,  in  the  sense  that  if  node  i  can  transmit  information  to 
another  node  j,  it  is  not  necessarily  true  that  node  j  can  transmit  information  to  node  i.  We  only  require 
that  each  node,  apart  from  seeing  incoming  transmissions  sent  to  it  by  neighboring  nodes,  knows  the 
number  of  nodes  that  it  can  transmit  information  to,  which  in  graph-theoretic  terms  is  referred  to  as  the 
out-degree  of  that  node.  In  fact,  in  the  proposed  algorithm,  each  node  will  broadcast  the  same  quantity  to 
all  receiving  nodes,  which  simplifies  the  communication  scheme  between  sending  and  receiving  nodes 
(as  it  is  not  necessary  for  each  sending  node  to  separately  communicate  with  each  receiving  node). 

When  the  communication  network  is  perfectly  reliable  (no  packet  drops),  the  collective  dynamics  of 
the  linear  iterations  can  be  described  by  a  discrete-time  transition  system  with  no  inputs  in  which  the 
transition  matrix  is  column  stochastic  and  primitive.  Then,  each  node  will  run  two  identical  copies  of 
the  linear  iteration  each  of  which,  however  is  initialized  differently  depending  on  the  problem  to  be 
solved.In  this  paper  we  mostly  focus  on  the  average  consensus  problem.  Under  proper  initialization, 

it  can  be  shown  that  each  node  will  asymptotically  calculate  the  desired  value  as  a  function  of  the 

outcomes  of  the  two  iterations.  The  details  of  these  double-iteration  approach  are  provided  in  [16], 

[17]  for  the  average  consensus  case  and  in  [8],  [9]  for  the  resource  coordination  problem.  For  the 

average-consensus  problem,  the  double-iteration  algorithm  is  a  particular  case  of  the  algorithm  in  [21] 
(which  is  a  generalization  of  the  algorithm  proposed  in  [22]),  where  the  matrices  describing  each  linear 
iteration  are  allowed  to  vary  as  time  evolves,  whereas  in  our  setup  (for  the  ideal  case  when  there  are 
no  communication  link  failures)  the  transition  matrix  is  fixed  over  time. 

The  focus  of  this  paper  is  to  robustify  the  double-iteration  algorithm  (informally  described  above  and 
formally  described  in  Section  II)  so  that  it  can  tolerate  failures  in  communication  links  and  converge 
to  the  average  value.  Our  communication  link  reliability  model  assumes  that  at  each  time  step,  a 
communication  link  is  unavailable  with  some  probability.  In  other  words,  a  packet  containing  information 
from  node  i  to  node  j  is  dropped  with  some  probability.  Next  we  informally  describe  our  robustification 
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approach.  Consider  two  nodes  i  and  j,  and  assume  that  j  receives  information  from  node  i  but  not 
necessarily  vice-versa.  Let  us  refer  to  j  as  the  receiving  node  (or  receiver)  and  i  as  the  sending  node 
(or  sender).  An  important  requirement  is  for  the  graph  describing  the  communication  network  to  be 
strongly  connected,  which  implies  every  node  must  be  able  to  act  both  as  a  sender  and  as  a  receiver. 
Then,  for  each  of  the  two  iterations  each  node  performs,  node  i  (the  sender)  will  keep  track  of  the 
following  quantities  of  interest:  i)  its  own  internal  state  (as  captured  by  the  state  variables  maintained  in 
the  double  iteration  scheme  of  [21],  [18];  ii)  the  total  mass  broadcasted  so  far  (to  be  described  in  detail 
soon);  and  iii)  the  total  received  mass  from  each  node  l  that  sends  information  to  node  i.  Similarly, 
for  both  iterations,  each  node  j  (the  receiver)  updates  the  value  of  its  internal  state  to  be  a  linear 
combination  of  its  own  previous  internal  state  value  (weighted  by  the  inverse  of  the  number  of  nodes 
that  have  j  as  a  neighbor)  and  the  difference  between  the  two  most  recently  received  mass  values  from 
each  of  its  neighbors  (also  weighted  by  the  inverse  of  the  number  of  nodes  that  have  j  as  a  neighbor). 
At  time  instant  k,  the  total  broadcasted  mass  by  node  j  is  the  sum  up  to  (and  including)  time  step  k  of 
the  sequence  of  values  of  node  j’s  internal  value,  weighted  by  the  inverse  of  the  number  of  nodes  that 
receive  values  from  node  j).  Additionally,  node  j  updates  the  value  of  the  received  mass  from  each  node 
l  that  sends  information  to  node  j  as  follows:  the  received  mass  from  node  l  is  the  total  broadcasted 
mass  sent  by  node  l  up  to  time  k  if  the  communication  link  from  node  l  to  node  j  is  available  at  time 
step  k;  otherwise,  the  received  mass  remains  the  same  as  the  most  recently  received  mass  from  node 
/.  An  implicit  assumption  here  is  that  messages  broadcasted  by  node  l  are  tagged  with  the  sender’s 
identity  so  that  the  receiving  node  j  can  determine  where  different  packages  have  originated  from. 

Recent  work  that  has  addressed  the  consensus  and  average-consensus  problems  in  the  presence  of 
unreliable  communication  links  [23],  [24],  [25]  has  employed  a  communication  link  availability  model 
similar  to  ours.  The  work  in  [23]  assumes  that  the  graph  describing  the  communication  network  is 
undirected  and  that  when  a  communication  link  fails  it  affects  communication  in  both  directions. 
Additionally,  nodes  have  some  mechanism  to  detect  link  unavailability  and  compensate  for  it  by  rescaling 
their  other  weights  (so  that  the  resulting  transition  matrix  remains  column  stochastic).  Following  this 
strategy,  the  authors  show  asymptotic  convergence  to  the  average  of  initial  conditions  and  also  calculate 
the  rate  at  which  the  variance  of  the  total  deviation  from  the  average  converges  to  zero.  The  work  in 
[24]  does  not  require  the  graph  describing  the  communication  network  to  be  undirected  and  proposes 
two  compensation  methods  to  account  for  communication  link  failures.  In  the  first  method,  the  so-called 
biased  compensation  method,  the  receiving  node  compensates  for  the  unavailability  of  an  incoming  link 
by  adding  the  weight  associated  to  the  unavailable  link  to  its  own  weight  (so  that  the  resulting  matrix 
remains  row  stochastic).  In  the  second  method,  called  the  balanced  compensation  method,  the  receiving 
node  compensates  for  link  unavailability  by  rescaling  all  the  incoming  link  weights  so  that  the  resulting 
matrix  remains  row  stochastic.  The  key  in  both  methods  is  the  fact  that  at  each  time  step,  the  resulting 
weight  matrix  is  row  stochastic;  the  authors  show  that  the  nodes  converge  almost  surely  to  the  same 
value,  but  this  value  is  not  necessarily  the  average  of  the  initial  conditions.  The  work  in  [25],  which  does 
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not  require  the  communication  graph  to  be  undirected,  proposes  a  correction  strategy  that  corrects  the 
errors  in  the  quantity  (state)  iteratively  calculated  by  each  node,  so  that  the  nodes  obtain  the  average  of 
their  initial  values.  This  correction  strategy  is  based  on  each  node  maintaining  some  auxiliary  variable 
that  accounts  for  the  amount  by  which  node  i  changes  its  state  due  to  the  updates  from  its  neighbors, 
i.e.,  the  nodes  that  can  send  information  to  node  i.  For  their  strategy  to  work  and  ensure  that  the  nodes 
converge  almost  surely  to  the  average  consensus,  the  authors  rely  on  the  nodes  sending  acknowledgment 
messages  and  retransmitting  information  an  appropriate  number  of  times. 

In  [21],  the  authors  proposed  a  gossip-based  algorithm  for  average-consensus  ver  a  directed  graph 
where  the  transition  matrices  describing  the  nodes’  collective  dynamics  change  at  every  iteration  step 
(depending  on  which  node  awakes).  This  scheme  requires  the  node  that  is  awake  to  perform  an 
internal  state  update  and  send  its  internal  state  (weighted  by  the  corresponding  out-going  link  weight) 
to  its  neighbors.  This  approach  results  in  generates  a  sequence  of  column  stochastic  matrices  (not 
necessarily  primitive)  with  the  property  that  all  the  diagonal  entries  remain  positive.  The  authors 
prove  that  by  running  two  such  iterations  in  parallel,  one  of  them  initialized  with  the  values  on 
which  the  average  operation  is  to  be  performed  and  the  other  with  the  all-ones  vector,  each  node 
will  asymptotically  achieve  average  consensus  by  taking  the  ratio  of  the  two  values  in  maintains.  A 
key  premise  in  their  proof  is  that  column  stochasticity  of  the  transition  matrix  is  maintained  over  time, 
which  requires  sending  nodes  to  know  the  number  of  nodes  that  are  listening.  This  suggests  that  i)  either 
the  communication  links  are  perfectly  reliable,  or  ii)  there  is  some  acknowledgment  and  retransmission 
mechanism  that  ensures  messages  are  delivered  to  the  listening  nodes  at  every  round  of  information 
exchange.  In  this  paper,  we  remove  such  assumptions  and  robustify  the  double-iteration  algorithm  against 
unreliable  communication  links  using  a  pure  broadcast-message  model  without  any  requirement  for  an 
acknowledgment/retransmission  mechanism.  Thus,  despite  the  reliance  of  our  algorithm  on  the  ratio  of 
two  linear  iterations,  it  is  different  both  in  the  communication  model  we  assume — a  broadcast  model 
in  our  case — and  also  in  the  nature  of  the  protocol  itself — our  focus  is  on  ensuring  convergence  in  the 
presence  of  communication  link  failures. 

An  additional  assumption  made  in  [21]  is  that  the  diagonal  entries  of  the  transition  matrix  (at  every 
step)  remain  positive.  In  our  model,  we  originally  consider  that  nodes  do  not  drop  self-packets.  However, 
to  ease  the  analysis,  we  remove  this  assumption  and  consider  the  case  where  self-packet  drops  are  also 
allowed  at  every  time  step,  which  i)  allows  us  to  handle  intermittent  faults  in  the  node  processing 
device,  and  ii)  removes  the  assumption  that  all  diagonal  entries  must  be  positive  at  every  step.  Finally, 
the  analysis  machinery  in  [21]  is  quite  different  from  the  one  used  in  this  paper.  We  employ  moment 
analysis  of  the  two  iterations  to  establish  that  they  are  linearly  related  as  the  number  of  steps  goes  to 
infinity,  while  [21]  relies  on  establishing  weak  ergodicity  of  the  product  of  the  transition  matrices  as 
the  number  of  steps  goes  to  infinity.  Finally,  as  it  will  be  shown  in  the  second  part  of  this  paper,  our 
algorithm  can  be  re-casted  into  a  similar  framework  as  the  one  in  [21]  by  augmenting  the  dimension 
of  the  vector  describing  the  collective  dynamics  to  account  for  the  packets  that  get  dropped  once  there 
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is  a  communication  failure.  Note,  However,  that  the  resulting  matrices  will  be  column  stochastic  but 
will  not  necessarily  have  strictly  positive  entries  on  their  diagonals.  In  the  second  part  of  the  paper,  we 
provide  an  analysis  framework  to  establish  the  convergence  of  our  algorithm  and  generalizes  the  ideas 
in  [21]  to  the  case  when  the  matrices  describing  the  system  collective  dynamics  do  not  have  strictly 
positive  diagonals.  In  this  regard,  we  will  show  that  even  in  the  case  where  self-packet  drops  are  not 
allowed,  the  resulting  transition  matrices  might  still  have  zero  diagonal  entries. 

The  remainder  of  this  paper  is  organized  as  follows.  Section  II  provides  background  on  graph  theory, 
introduces  the  communication  model,  and  briefly  describes  the  non-robust  version  of  the  double-iteration 
algorithm  we  use  in  this  work.  Section  III  describes  the  proposed  strategy  to  robustify  the  double-iteration 
algorithm  against  communication  link  failures  and  illustrates  the  use/performance  of  the  algorithms  via 
several  examples.  The  convergence  analysis  of  the  robustified  double-iteration  algorithm  is  provided  in 
Section  V.  Concluding  remarks  are  presented  in  Section  VI. 


II.  Preliminaries 


This  section  provides  background  of  graph-theoretic  notions  that  are  used  to  describe  the  communica¬ 
tion  network  and  the  distributed  consensus/coordination  setup,  introduces  the  basic  communication  link 
availability  model,  and  reviews  a  previously  proposed  two-iteration  algorithm  that  can  be  used  to  solve 
the  class  of  problems  addressed  in  this  paper  when  the  communication  network  is  perfectly  reliable. 

A.  Network  Communication  Model 

Let  discrete  time  instants  be  indexed  k  =  0, 1, ... ;  then,  the  information  exchange  between  nodes 
(components)  at  each  time  instant  k  can  be  described  by  a  directed  graph  Q[k]  =  (V,  £[k]},  where 
V  =  {l,2,...,n}  is  the  vertex  set  (each  vertex — or  node — corresponds  to  a  system  component),  and 
£[k\  C  V  x  V  is  the  set  of  edges,  where  (j,i)  G  £[k]  if  node  j  can  receive  information  from  node  i  at 
instant  k.  It  is  assumed  that  £  [k\  C  £,  \/k  >  0,  where  £  is  the  set  of  edges  that  describe  all  possibly 
available  communication  links  between  nodes;  furthermore,  the  graph  (V,  £)  is  assumed  to  be  strongly 
connected.  All  nodes  that  can  possibly  transmit  information  to  node  j  are  called  its  in-neighbors,  and 
are  represented  by  the  set  Af~  =  {i  G  V  :  (j,i)  G  £}.  Note  that  there  are  self-loops  for  all  nodes  in  Q 
(i.e.,  (j,j)  G  £  for  all  j  G  V).  The  number  of  neighbors  of  j  (including  itself)  is  called  the  in-degree 
of  j  and  is  denoted  by  Vj  =  \Nj~\-  The  nodes  that  have  j  as  neighbor  (including  itself)  are  called  its 
out-neighbors  and  are  denoted  by  AfJ~  =  {l  G  V  :  (/,  j)  G  £};  the  out-degree  of  node  j  is  P  ■  =  \Af+ 1 . 

The  existence  of  a  communication  link  from  node  i  to  node  j  can  be  described  in  probabilistic  terms 
as  follows.  At  instant  k,  let  Xji[k],  Vi,  j  G  V  be  an  indicator  variable  that  takes  value  1  with  probability 
q  and  takes  value  zero  with  probability  1  —  q,  i.e.. 


(1) 
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Then,  for  all  k  >  0,  the  existence  of  a  communication  link  between  nodes  i  and  j  can  be  described  be 
another  indicator  variable  £ji[k]  defined  as 


p  rtl  _  t  xAkl  if  CM)  e 

1  0,  if 

It  follows  that  £[k]  contains  the  elements  of  £  for  which  iyi [k]  =  Xji[k }  =  1. 


(2) 


B.  Double-Iteration  Algorithm  Formulation  in  Perfectly  Reliable  Communication  Networks 

When  the  communication  network  of  a  multi-component  system  is  perfectly  reliable,  i.e.,  Pr {£ji\k\  = 
1}  =  1,  G  £,  V/c  >  0,  it  was  shown  in  [17],  [9]  that  the  components  of  the  multi-component 

system  can  asymptotically  solve  average  consensus  and  resource  coordination  problems  in  a  distributed 
fashion  by  running  two  separate  appropriately  initialized  linear  iterations  of  the  form 

vAk  +  !]  = 

zAk  +  !]  =  5^  ^+zAkl  (4) 

i&A/T  4 


where  T>k  (Vf  )  is  the  out-degree  of  node  j  (/’).  A  requirement  in  all  cases  is  that  the  underlying 
communication  graph  (Q,£)  is  strongly  connected. 

1 )  Average  Consensus  Problem:  In  this  problem,  the  nodes  aim  to  obtain  the  average  of  the  values 
Vj ,  j  =  1, ...  ,n,  they  each  posses.  In  [17],  it  was  shown  that  if  the  initial  conditions  in  (3)  (referred 
to  as  iteration  1)  are  set  to  yf  0]  =  vv  and  the  initial  conditions  in  (4)  (referred  to  as  iteration  2)  are 
set  to  Zj[ 0]  =  1,  then  the  nodes  can  asymptotically  calculate  v  :=  J2'j=ivj/n  as 


v  = 


lim 

k — ^oo 


vAk ] 
zAkY 


(5) 


by  running  the  two  iterations  in  (3)  and  (4). 

2)  Resource  Coordination  Problem:  In  this  problem,  each  node  j  can  contribute  a  certain  amount 
7 Tj  >  0  of  a  given  resource,  which  is  upper  and  lower  bounded  by  known  capacity  limits  N™x  and 
7r"nn  respectively.  The  challenge  is  to  coordinate  the  components  so  that  they  collectively  provide  a 
pre-determined  total  amount  pd  =  V"  ,  7 Tj  of  the  resource1  as  specified  by  an  external  “leader.”  In  [9], 
it  was  shown  that  i)  if  the  initial  conditions  in  (3)  are  set  to  yA^\  —  f>d/rn  ~  n™7'  if  j  is  an  out-neighbor 
of  the  leader  (where  m  >  1  is  the  number  of  nodes  contacted  initially  by  the  external  leader)  and 
Dj[ 0]  =  —  7r’nm  otherwise,  and  ii)  if  the  initial  conditions  in  (4)  are  set  to  Zj [0]  =  —  tt""".  then  the 


'in  the  development  in  [18],  [9],  it  is  assumed  that  XP=i  <  pd.  <  57"=i  ’K1Jla'x\  this  is  not  a  restrictive  assumption  because  in  the 
proposed  algorithms,  all  nodes  will  be  able  to  know  if  pd  <  f£j=i  or  if  Pd  >  57^=1  nj1’ax  (which  means  that  no  feasible  solution 
exists). 
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nodes  can  asymptotically  calculate  their  own  resource  contribution  7 Tj  as 

yAk\ 


Ti :=  lim  ( 

k—>  00 


rnin  .  L  J  /max  rnin 

+W)(i  ~ 1 


\\  _  m 

) )  — 77 j 


En  r 

^  1=1 77 1 

\  . . "  l 


Eiu*. 


/  ^.max  rnin 


), 


(6) 


which  satisfies 


rnin  ^  ^  max  \_v  • 

77  j  <  77  j  < 77  j  ,  Vj, 

n 

=  Pd •  (7) 

7=1 

In  this  paper,  we  start  with  a  double  iteration  of  the  form  in  (3)-(4)  that  is  used  for  either  aver¬ 
age  consensus  or  coordination,  and  develop  systematic  methodologies  to  handle  packet  drops  in  the 
communication  links. 

III.  Robustification  of  Double-Iteration  Algorithm 

In  this  section,  the  algorithm  described  in  Section  II-B  is  modified  so  as  to  make  it  robust  against 
communication  link  failures.  As  in  (3)-(4),  each  node  will  run  two  iterations  (which  we  refer  to  as 
iterations  1  and  2)  to  calculate  quantities  of  interest  and  eventually  solve  the  average  consensus  or 
coordination  problems. 

Consider  the  setup  described  in  the  previous  section:  we  are  given  a  (possible  directed)  strongly 
connected  graph  (Q .  £)  representing  a  multicomponent  system  and  its  communication  links  between  its 
components.  For  the  sake  of  generality,  let  us  refer  to  j  as  the  receiving  node  (or  receiver)  and  i  as 
the  sending  node  (or  sender).  For  each  of  the  two  iterations,  node  i  (the  sender)  will  calculate  several 
quantities  of  interest,  which  we  refer  to  as:  i)  internal  state;  ii)  total  broadcasted  mass;  and  iii)  total 
received  mass  from  each  in-neighbor  l  of  node  i,  i.e.,  for  each  node  l  G  .  For  both  iterations  1  and  2, 
each  node  j  updates  the  value  of  its  internal  state  to  be  a  linear  combination  of  its  own  previous  internal 
state  value  (weighted  by  the  inverse  of  the  number  of  nodes  that  have  j  as  a  neighbor,  i.e.,  l/Pt)  and 
the  sum  (over  all  its  in-neighbors)  of  the  difference  between  the  two  most  recently  received  mass  values. 
At  instant  time  k,  the  total  broadcasted  mass  is  the  sum  up  to  (and  including)  step  k  of  the  weighted 
value  of  node  j’s  internal  state  (used  to  update  the  internal  state  of  node  j).  Additionally,  node  j  (the 
receiver)  updates  the  value  of  the  received  mass  from  node  /  to  be  either  the  total  broadcasted  mass  sent 
by  node  i  if  the  communication  link  from  i  to  j  is  available  at  instant  k,  or  the  most  recently  received 
mass  value  from  node  i,  otherwise.  An  implicit  assumption  here  is  that  messages  broadcasted  by  node 
i  are  tagged  with  the  sender’s  identity  so  that  the  receiving  node  j  can  determine  where  messages 
originated  from. 

For  iteration  1,  let  yj[k]  denote  node  j' s  internal  state  at  time  instant  k,  jiij [/;']  denote  the  mass 
broadcasted  from  node  j  to  each  of  its  out-neighbors  l  (this  is  a  single  value  ad  it  is  the  quantity  is 
the  same  for  each  out-neighbor  l  of  node  j,  i.e.,  for  each  l  e  A^A,  and  vjr  [/,: |  denote  the  total  mass 
received  at  node  j  from  node  i  E  J\f~ .  Similarly,  let  z3  \k]  denote  node  j’s  internal  state  takes  at  time 
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instant  k,  crij[k\  denote  node  j' s  broadcasted  mass  for  each  out-neighbor  l,  l  E  N-  ,  and  t3% [k]  denote 
the  total  mass  received  from  node  i  E  .  Then,  the  progress  of  iteration  1  is  described  by 

vAk  + !]  =  ^+vAk)  +  feM  -  uAk  -  !])>  k  > 

3  i&A/T 

1  '  ]_ 

imAA  =  mAk  - !]  +  =  J]  j^vA'A  k  >  o,  (8) 

i= 0  U3 


where 


VjAk\ 


fiji[k\,  if  (j,i)  E  S[k],  k>  0, 

Uji[k  -  1],  if  (J,i)  £  £[k\,  k>  0. 


Recall  that  V t  (X>+)  is  the  number  of  nodes  that  node  j  (i)  can  transmit  information  to.  Similarly,  the 
progress  of  iteration  2  is  described  by 


ie  M. 


Vi 


[vAk\  -  vAk  - !]), 

k  >  0, 

k 

#]  = 

i= 0  ui 

k  >  0, 

(9) 

where 


Tji[k\  = 


CTji[k\ ,  if  (j,i)  E  £[k],  k>  0, 

Tji[k  -  1],  if  (j,i)  £  £[k],  k>  0. 


As  mentioned  earlier,  for  solving  the  average  consensus  problem,  the  initial  conditions  in  (8)  are 
set  to  yj[ 0]  =  Vj,  whereas  the  initial  conditions  in  (9)  are  set  to  z3 [0]  =  1.  Similarly,  for  solving  the 
resource  coordination  problem,  the  initial  conditions  in  (8)  are  set  to  y3  [0]  =  pci /rri  —  7iJan  if  j  is  initially 
contacted  by  the  leader  and  y3  [0]  =  —  tt"""  otherwise,  whereas  the  initial  conditions  in  (9)  are  set  to 
Zj[ 0]  =  7iJMX  —  it ™-in  >  0.  In  both  the  average  consensus  and  coordination  problems,  //T,[— 1]  =  0  and 
Vji[—  1]  =  0  for  all  (J,i)  E  £,  and  <J3i[—l]  =  0  and  r3i[—l]  =  0  for  all  (j,  i)  E  £. 

Main  Result:  We  shall  argue  that  with  the  proposed  robustification  strategy,  despite  the  presence  of 
unreliable  communication  links  (at  each  time  step,  each  link  (j,  i )  E  £,  fails  independently  from  other 
links  and  independently  between  time  steps,  with  some  probability  1  —  q3l),  nodes  can  asymptotically 
estimate  the  exact  solution  v  to  the  average  consensus  by  calculating,  whenever  z3  [k]  >  0  the  ratio 
yAk\/zAk\,  i.e., 


v  = 


lim 

k—>  oo 


vAk\ 

ZjW 


(10) 
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whenever  Zj [k  >  0|.  Similarly,  exact  solution  to  the  resource  coordination  problem  can  be  obtained  as 

_  ( ^.min  ,  Ujl^i  ( ^.max  _ min \\  /i  i\ 

+wi(  i  1  )]' 

whenever  Zj[k\  >  0.  In  both  cases,  we  run  the  iterations  in  (8)  and  (9)  and  using  the  corresponding  initial 

y^n_  zo(j) 

conditions  as  outlined  above.  In  particular,  we  will  show  that,  for  every  j,  zj[k\  —  jjp 1  y^k) y3 [k]  — >  0 
as  k  — *  oo  almost  surely.  Additionally,  we  will  show  that  z3  [k]  >  0  occurs  infinitely  often. 

A.  Examples 

We  now  illustrate  how  the  proposed  algorithm  works  for  the  case  of  average  consensus  in  the 
presence  of  packet-dropping  communication  links.  We  start  with  the  rather  small  network  shown  in 
Fig.  1  and  assume  that  the  packets  on  each  link  (including  the  self-links  which  are  not  drawn  in  the 
figure2)  can  be  dropped  with  probability  1  —  q,  independently  between  different  links  and  independently 
between  different  iterations.  We  also  assume  that  the  initial  values  of  the  five  nodes  are  given  by 
v  =  [—4,  5,  6,  —3, 1]T,  with  their  average  equal  to  1.  Thus,  in  the  iterations  (8)  and  (9) 

y[ 0]  =  [-4,  5,  6,  -3, 1]T,  and  z[0]  =  [1, 1, 1, 1, 1]T  , 

with  pji[- 1]  =  Vji[-1]  =  (Tji[- 1]  =  Tji[- 1]  =  0  for  all  (j,i)  e  E. 

We  run  the  iterations  in  (8)  and  (9)  and  plot  the  ratio  |4^|  as  a  function  of  the  iteration  step  k  for  each 
node  j  ( j  =  1,  2,  3, 4,  5).  Figure  2  shows  the  typical  behavior  that  we  observe  for  q  =  0.99  (i.e.,  for  a 
probability  of  a  packet  drop  equal  to  0.01).  As  can  be  seen  in  the  figure,  the  ratio  at  each  node  quickly 
converges  to  the  correct  average,  though  the  individual  values  for  y3  [k]  and  z3  [k]  do  not  converge. 

In  Figs.  3  and  4  we  show  typical  behaviors  of  the  same  multicomponent  system  for  q  =  0.5  and 
q  =  0.1.  The  behavior  remains  similar  to  the  one  observed  before:  even  though  y3 [k]  and  z3  [k]  do  not 

2We  make  this  assumption  later  in  the  paper  for  the  purposes  of  simplifying  notation. 


Fig.  1.  Small  directed  graph  used  for  illustration  of  the  ratio  algorithm  for  obtaining  average  consensus  in  the  presence  of  packet-dropping 
communication  links. 
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y  at  Each  Node  vs  Iterations 


z  at  Each  Node  vs  Iterations 


Number  of  Iterations 


Ratio  at  Each  Node  vs  Iterations 


Number  of  Iterations 


(a) 


(b) 


Fig.  2.  Evolution  of  the  values  of  yj[k\  (left),  Zj[k\  (middle)  and  l^jjq  (right)  for  q  =  0.99. 


(c) 


Ratio  at  Each  Node  vs  Iterations 


Number  of  Iterations 


(a) 


(b) 


(c) 


Fig.  3.  Evolution  of  the  values  of  j /j  [fc]  (left),  Zj  [fc]  (middle),  and  ^-||j  (right)  for  q  =  0.5,  j  =  1,2,3, 4,  5. 


converge  (in  fact,  they  seem  to  behave  more  radically  with  decreasing  q),  the  ratio  |4^|  does  converge 
to  the  average  of  the  initial  values.  Note  that  the  plot  of  the  ratio  in  Fig.  4  is  quite  different  than  the 
rest:  in  this  case,  q  is  small  enough  so  that  Zj[k]  (and  simultaneously  yj[k])  can  take  the  value  zero 
(e.g.,  when  all  packets  destined  to  node  j  are  dropped  at  iteration  k );  thus,  the  ratio  in  such  cases  is 
not  defined  and  is  not  plotted,  resulting  in  a  discontinuous  set  of  points  in  the  plot.  Nevertheless,  we 
can  see  that  when  packets  are  received  (which  happens  frequently  enough  for  each  node),  the  ratio  has 
the  correct  value.  This  is  a  point  addressed  later  in  the  paper. 

An  example  of  what  happens  in  larger  graphs  is  shown  in  Fig.  5.  Here  we  consider  a  graph  with  50 
nodes,  randomly  generated  by  choosing  a  directed  edge  from  node  i  to  node  j ,  1  <  i,j  <  50,  i  ^  j , 
independently  with  probability  1/2,  and  ensuring  that  the  resulting  graph  is  strongly  connected.  As  can 
be  seen  the  behavior  remains  similar  to  what  we  observed  for  the  smaller  graph:  the  ratio  converges 
quickly  to  the  average  even  though  the  individual  y3 [k]  and  Zj[k\  do  not  converge.  For  this  particular 
plot,  we  used  q  —  0.1,  which  also  justifies  the  fluctuation  in  the  values  of  y,  [k]  and  z3  [k] . 
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Ratio  at  Each  Node  vs  Iterations 


(a) 


(b) 


(c) 


Fig.  5.  Evolution  of  the  values  of  yj[k]  (left),  Zj[k\  (middle),  and  ^-||j  (right)  for  q  =  0.1  for  a  50-component  system. 


IV.  First  and  Second  Moment  Analysis 

In  this  section,  we  obtain  recurrence  relations  that  describe  the  first  and  second  moment  of  the 
iterations  after  (8)  and  (9);  this  analysis  is  used  in  Section  V  to  establish  the  claims  in  (10)  and  (11).  In 
order  to  ease  the  moment  calculations,  the  expressions  in  (8)— (9)  will  be  rewritten  more  compactly  in 
vector  form.  Also,  in  order  to  facilitate  notation,  we  will  allow  each  node  j  to  drop  the  packet  carrying 
its  own  previous  value  when  updating  its  value.  This  way,  node  j  handles  its  own  value  in  the  same 
way  as  its  neighbors’  values  and  notation  is  simplified  significantly. 

A.  Vectorized  Description  of  Double-Iteration  Algorithm 

Using  the  definition  for  the  indicator  variable  x3,\k\  given  in  (1)  and  the  resulting  indicator  variable 
£ji[k\  given  in  (2),  which  describes  the  successful  transmission  of  information  from  node  i  to  node  j 


COORDINATED  SCIENCES  LABORATORY  TECHNICAL  REPORT  UILU-ENG-1 1-2207  (CRHC- 11-05) 


13 


over  an  existing,  unreliable  communication  link,  iterations  (8)  and  (9)  can  be  rewritten  as 


and 


Vjilk] 


Mi [k  ~  !]  +  Vj[k],  iflett’f,  k>  0, 

0,  if  l<£  A f+,  k  >  0, 

liji[k\xji[k\  +  Uji[k  -  1](1  -  Xji[k]),  if  i  G  Af~,  k>  0, 

0,  if  i  ^  A f~,  k  >  0, 


Vj[k  +  1]  =  (; vji[k }  -  Uji[k  -  1]),  k>  0, 

i=l 


aij  [k] 


Tji[k} 


crij[k  -  1]  +  ^+Zj[k\,  if  /  G  Af+,  k  >  0, 

0,  if  Z  ^  Aff,  k>  0, 

aji[k]xji[k\  +  Tji[k  -  1] (1  -  Xji[k}),  if  i  G  Af~ ,  k  >  0, 

0,  if  i  G  Af~,  k  >  0, 


Zj[k  +  1]  =  ( Tji[k\  -  Tji[k  -  1]),  k>  0, 

l— 1 


(12) 

(13) 

(14) 

(15) 

(16) 

(17) 


where  ptj[- 1]  =  =  <7y[-l]  =  ^i[— !]  =  0,  Vj,i. 

Let  Ao  B  denote  the  Hadamard  (entry-wise)  product  of  a  pair  of  matrices  A  and  B  of  identical  size. 
Then,  for  all  k  >  0,  iteration  ( 1 2)— ( 14)  can  be  rewritten  in  matrix  form  as 


Mk  =  Mfc_!  +  Pdiag(yfc),  (18) 

Nk  =  MkoXk  +  Nk_ !  o  (U  -  Xk),  (19) 

Vk+i  =  (Nk  ~  Nk-i)e  =  \(Mk  -  Nk _x)  o  Xk\  e,  (20) 

where  P  =  \pji ]  G  Mnxn,  with  pji  =  ,  Vj  G  A ft+  and  pyi  =  0  otherwise;  M_i  =  iV_ i  =  0 \yk  =  y[k\, 
U  G  M”'xn,  with  [Uji\  =  1,  Vi,  j;  diag(y/.)  is  the  diagonal  matrix  that  results  by  having  the  entries  of 
yk  on  the  main  diagonal;  and  e  =  [1, 1, . . . ,  l]7  (note  that  U  =  eeT).  Similarly,  for  k  >  0,  (15)— (17)  can 
be  rewritten  in  matrix  form  as 


Sk  =  Sk~i  +  Pdiag(^),  (21) 

Tk  =  Sk  o  Xk  +  Tfc_ i  o  (U  —  Xk),  (22) 

zk+ 1  =  (Tk  —  Tk_ i)e  =  [(Sk  —  Tk-i)  o  Xk]e,  (23) 

where  5_ i  =  T  \  =  0,  =  z [/;•] ,  and  diag(zfc)  is  the  diagonal  matrix  that  results  by  having  the  entries 

of  zk  on  the  main  diagonal. 

By  defining  Ak  :=  Mk  —  Nk_\  and  Bk  :=  Sk  —  Tfe_i,  iteration  (18)-(20)  can  be  rewritten  more 
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compactly  as 


Ak  =  o  (U  +  Pdiag(yk),  k  >  1,  (24) 

Z/fc+i  =  (Ak  o  Xfc)e,  k  >  0,  (25) 

and  iteration  (21)— (23)  as 

Bk  =  Pfe_i  o  (U  -  Xfe_i)  +  Pdiag(zfe),  k  >  1,  (26) 

£fc+i  =  (-Bfc  o  Xfe)e,  k  >  0,  (27) 

where  /10  =  M0  —  N_}  =  Pdiag(i/0),  and  B0  =  S0  —  T_i  =  Pdiag(z0). 

For  analysis  purposes,  each  matrix  in  (24)-(25)  and  (26)-(27)  will  be  rewritten  in  vector  form  by 
stacking  up  the  corresponding  columns.3  Then,  (24)-(25)  and  (26)-(27)  can  be  rewritten  in  vector 
form  as  follows.  Let  F  =  [In  In  ...  In\  E  Rnxn2,  where  In  is  the  n  x  n  identity  matrix,  and  P  = 
[ EiPT  E2Pt  . . .  EnPT]T  G  Wi"2xn,  where  E,  e  Mnx"  has  E,(i-  i)  =  1  and  all  other  entries  equal  zero. 
[The  entries  of  EiPT  e  Mnx"  (PEf  =  PEi)  are  all  zero  except  for  the  ith  row  (column)  entries,  which 
are  those  of  the  ith  row  (column)  of  matrix  PT  (P).]  Then,  (24)-(25)  can  be  rewritten  as 

ak  =  ak_ i  o  (u  -  xk-i)  +  Pyk,  k>  1,  (28) 

Vk+i  =  F(ak  o  xk),  k>  0,  (29) 

2  2  2 

where  ak  E  Mn  ,  xk  £  Mn  ,  and  xk-\  E  Mn  result  from  stacking  the  columns  of  matrices  Ak,  Xk,  and 
Xk_i,  respectively.  Similarly,  (26)-(27)  can  be  rewritten  as 

bk  =  bk- !  o  (u  -  xk-i)  +  Pzk,  k>  1,  (30) 

zk+i  =  F(bk  o  xk),  k>  0,  (31) 

where  bk  E  M”  results  from  stacking  the  columns  of  matrix  Bk. 

Remark  1:  It  is  important  to  note  that  matrices  Ak  and  Bk,  and  their  corresponding  vectors  ak  and 
bk,  have  some  entries  that  remain  at  zero  for  all  k  >  0.  Specifically,  the  (),  i)  entry  of  matrices  Ak  and 
Bk  (and  their  corresponding  entries  in  ak  and  bk)  remain  zero  if  there  is  no  communication  link  from 
node  i  to  node  j,  i.e.,  (j,  i )  ^  £.  The  reason  we  keep  these  entries  (despite  the  fact  they  are  zero  and  do 
not  play  a  role  in  the  analysis)  is  because  they  facilitate  matrix  notation  and  calculations  in  subsequent 
developments.  □ 

Since  it  will  appear  later  at  several  points  of  the  analysis,  it  is  worth  noting  that  when  premultiplying 
P  by  F,  we  recover  the  matrix  P,  i.e., 


P  =  PP. 


(32) 


’If  we  let  A  =  [ Aij ]  €  Rnxn,  then  a  =  [An,  A2i,  . . . ,  Am,  A12,  A22,  ■  ■  ■ ,  An2, . . . ,  Ain,  A2n,  ■  ■  ■  ,  Ann]T .  Vectors  defined  by 
stacking  the  columns  of  a  matrix  will  be  denoted  with  the  same  small  letter  as  the  capital  letter  of  the  corresponding  matrix. 
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B.  First  Moment  Analysis 

In  this  section,  we  describe  the  first  moment  dynamics  of  (28)— (3 1)  via  discrete-time  transition  systems 
with  no  inputs,  where  (as  shown  below)  the  corresponding  transition  matrices  (which  are  obtained  from 
P  and  q)  are  column  stochastic  and  primitive.  In  both  iterations,  the  sum  of  the  entries  of  the  first 
moment  vectors  for  yk  and  zk  is  shown  to  remain  constant  over  time  and  be  respectively  equal  to  the 
sum  of  qYliV  o(*)  and  q  zoH)-  Furthermore,  both  first  moments  E  [yk]  and  E[~fc]  are  shown  to  reach 
a  steady-state  value  as  k  goes  to  infinity.  The  above  discussion  is  formalized  in  the  following  lemma. 

Lemma  1:  Let  ak,  bk,  IJk,  and  Zk  be  described  by  the  recurrence  relations  in  (28)-(29),  and  (30)-(31) 
respectively.  Let  the  first  moments  of  ak,  Vk ,  bk,  and  Zk  (i.e.,  E[afc],  E [yk\,  E [bk],  and  E[^-])  be  denoted 
by  ak,  Vk,  bk,  and  zk  respectively.  Then  the  evolution  of  ak,  yk,  bk,  and  zk,  \/k  >  1,  is  governed  by 


ak  =  [qPF  +  (1  -  q)In2]ak-i, 

(33) 

yk+i  =  [qP  +  (i  -  q)in\yk, 

(34) 

bk  =  [qPF  +  (1  -  q)In2]bk-i, 

(35) 

zk+i  =  [qP  +  (1  -  q)In]zk, 

(36) 

where  Im  is  the  m  x  m  identity  matrix,  with  a0  =  Pyo,  if  =  qPyo,  b0  =  Pz0,  and  z\  =  qPz0. 

Proof:  Since  the  development  for  obtaining  ak  and  yk  is  parallel  to  that  for  obtaining  bk  and  zk, 
our  analysis  focuses  on  the  first  case.  For  k  =  0  in  (28)-(29),  by  taking  expectations  of  both  sides  and 
noting  that  packet  drops  at  time  step  k  —  0  are  independent  of  the  initial  values  for  ao,  it  follows  that 

a0  =  Pyo,  (37) 

y,  =  qFa0.  (38) 

Substituting  (37)  into  (38),  we  obtain  y1  =  qFPy0  =  qPy0- 
For  k  >  1  in  (28)-(29),  noting  that  packet  drops  at  time  step  k  are  independent  of  previous  packet 


drops  and  the  initial  values  of  a0,  it  follows,  by  taking  expectations  on  both  sides,  that 

ak  =  ak- 1  o  (u-  xk-i)  +  Pyk  =  ak- 1  o(u-  xk-i)  +  Pyk  =  (1  -  q)ak- 1  +  Pyk,  (39) 
yk+ i  =  F(ak  o  xk)  =  F(ak  o  xk)  =  qFak.  (40) 

Substituting  (40)  into  (39),  we  obtain 

ak  =  (1  -  q)ak~i  +  qPFcik- 1  (41) 

=  [ qPF  +  (1  -  q)In2\ak-\,  (42) 
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Similarly,  substituting  (39)  into  (40),  we  have 


Uk+ 1 

=  (1  -  q)qFak_i  +  qFPyk 

(43) 

=  (i  -  q)yk  +  q.FPyk 

(44) 

=  [qP  +  (!  -  d)in]yki 

(45) 

where  In  is  the  n  x  n  identity  matrix. 

■ 

C.  Second  Moment  Analysis 

In  the  order  to  calculate  the  second  moment  dynamics  for  (28)— (3 1),  we  utilize  in  the  following 
lemma. 

Lemma  2:  Let  x,  c  and  d  be  random  vectors  of  dimension  n.  Furthermore,  assume  that  the  entries  of  x 


are  Bernoulli  i.i.d.  random  variables  such  that  Pr{xj  =  1}  =  q  and  Pr{xj  =  0}  =  1  —  q,  Vi  =  1,  2, . . .  n, 
and  are  independent  from  c  and  d.  Then 

S  :=  E  [(c  o  x)(x  o  d)7}  =  q2  E [cdT]  +  q(  1  —  q)  E  [diag(ccF)],  (46) 

T  :=  E  [(co  x)  ((m  —  x)  o  d)T]  =  q(  1  —  q)  E [cdT]  —  q(  1  —  q)  E  [diag(cdT)],  (47) 

where  diag(ccF)  is  a  diagonal  matrix  with  the  same  diagonal  as  matrix  cdT. 

Proof:  The  (■ i,j),i  f  j,  entry  of  S  can  be  obtained  as  follows: 

Sij  =  E  [axidjXj] .  (48) 

Since  X;  and  Xj  are  pairwise  independent,  and  independent  from  c  and  d,  it  follows  that 

E  [axidjXj]  =  q-  E  [c,d:j\ .  (49) 


For  i  =  j,  observing  that  E [x;Xj]  =  E[x;]  =  q,  Vi  =  1, . . . ,  n,  we  obtain  the  corresponding  entry  of  S 
as 

Su  =  E  [c^Axi]  =  E  [cidiXi]  =  q  E  [cA] .  (50) 

In  (49),  it  is  easy  to  see  that  E  \cidj\  is  the  (i,j)  entry  of  E [cdT].  Similarly,  in  (50),  it  is  easy  to  see 
that  E  [cjdj]  is  the  (?',  i)  entry  of  E [cdT].  From  these  observations,  the  result  in  (46)  follows. 

Similarly,  the  ( i,j),i  f  j ,  entry  of  T  can  be  obtained  as  follows: 

Tij  =  E  [dXidjil  -  Xj )] .  (51) 

Since  X;  and  (1  —  xf)  are  independent,  it  follows  that 

E  [ciXidj(  1  -  Xj)]  =  E  [cidj]  E  [xj(l  -  Xj)]  =  q{  1  -  q)  E  [cidj] . 


(52) 


COORDINATED  SCIENCES  LABORATORY  TECHNICAL  REPORT  UILU-ENG-1 1-2207  (CRHC- 11-05) 


17 


For  i  =  j,  and  observing  that  E[x*(l  —  x*)]  =  0,  Vi  =  1, . . .  ,n,  the  corresponding  entry  of  T  can 
obtained  as  follows; 


Ta  =  E  [ciXidi{  1  -  xi)]  =  E  [ddi]  E  [xj(l  -  x*)]  =  0. 


(53) 


The  result  in  (47)  follows  from  (52)  and  (53).  ■ 

The  following  lemma  establishes  that  the  evolution  of  E [akak\,  E[bkb^\,  and  E[akbk],  and  can  be 
expressed  as  linear  iterations  with  identical  dynamics  but  different  initial  conditions.  Similarly,  the 
evolution  of  E[yky^],  E[zkzk],  and  E[ykzk]  can  also  be  expressed  as  linear  iterations  with  identical 
dynamics  but  different  initial  conditions. 

Lemma  3:  Consider  the  second  moments  of  ak,  yk,  bk,  and  zk,  and  let  E [akap\,  E [ykyk] ,  E [bkbk], 
Ej zkz[\,  E [akbp\,  and  E [ykzk])  be  denoted  by  Tk,  $fc,  T/,:,  A/,,  Ek,  and  T k  respectively.  Then,  the 
evolutions  of  Tk,  T/,:,  Afe,  Ek,  Tk.  Vk  >  1,  are  described  by  the  following  iterations  (where  all  / 
denote  n2  x  n2  identity  matrices): 

Tfc  =  [qPF  +  (1  -  g)J] rfc_!  [qPF  +  (1  -  q)l]T  +  q(  1  -  q)[I  -  P F]dmg{T ^[I  -  PF]T,  (54) 
<f>fc+i  =  F  [q2T k  +  q(  1  -  g)diag(T fe)]  FT,  (55) 

T,  =  [qPF  +  (1  -  q)l]  [qPF  +  (1  -  q)l] T  +  q(l  -  q)  [I  -  PF] diag(^fc_1)  [/  -  PF}T,  (56) 
Afc+i  =  F[q2' +  q(  1  -  g)diag('Tfc)]FT,  (57) 

Ek  =  [qPF  +  (1  -  g)/]-^!  [qPF  +  (1  -  q)l]T  +  q{  1  -  q)[I  -  PE]diag(Sfe_1)[J  -  PF]T,  (58) 
T k+i  =  F  [q2Ek  +  q(  1  -  g)diag(Sfc)]  PT,  (59) 

with  initial  conditions 


r0  =  Py0yoPT,  (60) 

=  y\VT\  +  ?(1  -  q)Fdiag(Py0yTpT)FT,  (61) 

'&o  =  Pzqz£Pt,  (62) 

Ai  =  zxz\  +  g(l  -  g)Fdiag(P^o^PT)FT,  (63) 

E0  =  Py0ZoPT,  (64) 

Ti  =  y{z\  +  g(l  -  q)Fdiag(Py0ZQ  PT)FT .  (65) 


Proof:  The  derivation  of  (54),  (55),  (60),  and  (61),  is  the  same  as  the  derivation  of  (56),  (57),  (62), 
and  (63),  thus  the  developments  in  the  proof  will  only  address  the  former.  For  k  —  0,  it  follows  from 
Lemma  1  and  (28)  that  ao  =  Pyo ■  Then, 


T0  =  E[a0ag  ]  =  PE[y0y^]PT  =  Py0y^  PT , 


(66) 
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and 


=  ElyxyJ]  =  E  [F(a0  o  x0)(x0  o  a0)TF7]  =  F  E  [(a0  o  x0)(x0  °  a0)T]F 1 .  (67) 


Applying  the  results  in  Lemmas  1  and  2  to  (67),  it  follows  that 

<f>i  =  q2F E  [aQal]FT  +  q{  1  -  q)F  E  [diag(a0ao  )j  FT 

=  ( qFPyo)(qFPy0)T  +  q(  1  -  q)F  E[diag(Py0y^  PT)\FT 

=  ( qPy0)(qPyo)T  +  q(  1  -  g)Fdiag(P7/o^PT)FT 

=  ViVl  +  7(1  -  7)^diag(Py0|/(rPr)PT,  (68) 

where  we  used  the  fact  that  FP  =  P  (refer  to  Eq.  (32)) . 

For  k  >  1,  and  taking  into  account  that  yk  =  F(ak_i  °  xk_ i),  it  follows  that 

=  E  [(ak-i  o(u—  xk-i)  +  Pyk)  (ak-i  °  (u  -  xfc_i)  +  Pyk)T] 

=  E  [(ak-i  o  (u-  xk-i))  (ak-i  o  (u  -  xfc_i))T]  +  E  [(afc_i  °  (it  Zfc-i))  (. Pyk)T ] 

+  E  [(Pyk)  (ajfc-i  o  (it  -  £fc-i))T]  +  E  [(Pyfc)  (Pyk)T] 

—  E  [(afc_i  o  (it  -xfc_i))(afc_i  o  (u-xk_i))T]  +E  [(afc_ i  o  (it  -  xfc_i))  (afc_i  o  xfc_i)r]  PTPT 

+  PPE  [(afc_i  oxfc_i)(afc_i  o  (u  -  xfc_i))T]  +  PPE  [(afc_i  o  xfc_i)  (afc_i  o  xk_i)T]  FTPT. 

(69) 

Then,  from  Lemma  2,  (69)  can  be  rewritten  as 

Efc  =(1  -  q )2  E  [afc_iafc_J  +  g(l  -  q)  E  [diag^^a^)] 

+  7(1  -  7)  E  [ak-1al_1]FTPT  -  q(  1  —  g)  E  [diag(afc_ia^_1)]PrPT 
+  7(1  -  7)P P  E  [ajfc_i a^_J  -  g(l  -  g)PFE  [diag(afc_ia^’_1)] 

+  g2P P  E  [afc_1a^’_1]PTPT  +  g(l  -  q)PF  E  [diag(afc_1a^_1)]FTPr.  (70) 

By  re-arranging  terms  in  (70)  and  observing  that  rfc_x  =  E  \ak_ia^_^  and  diag(rfc_!)  =  E  [diag(afc_ia|L1)] , 
the  result  in  (54)  follows. 

Additionally,  from  Lemma  2,  it  follows  that 


®k+i  =  E  [yk+iyl+i] 


F  E  [(ak  o  xk)  ( xk  o  ak)T]  F 7 


F 


q2  E  [akal]  +  q(  1  -q)E  [diag(afcaf )] 


F  [g2r  k  +  q(l  —  g)diag(r  fc)]  F7 


(71) 
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To  obtain  the  iterations  for  Ek  and  T k,  the  developments  are  very  similar  to  the  ones  above.  For 
k  —  0, 

50  =  E  [a0&o  ]  =  EE[t/o4]Pt  =  PVoZqP,  (72) 

and 


Ti  =  E[yizf]  =  E  [F(a0ox0)(x0  °  b0)rFT]  =  F E  [(a0oi0)(a;o  o  b0)T]FT 
=  (qFPy0)(qFPz0)T  +  q{  1  -  g)FE[diag(P^0TPT)]FT 
=  (qPyo){qPz0)T  +  q{  1  -  g)Fdiag(P7/0z[[PT)FT 

=  VFi  +  q{  1  -  q)Fdiag(Py0ZQpT)FT,  (73) 

where  again  we  used  the  fact  that  FP  =  P  (refer  to  Eq.  (32)) . 

For  k  >  1,  from  Lemma  2  and  (28),  and  taking  into  account  that  yk+i  =  F(ak  °  xk)  and  zk+ 1  = 
F(bk  o  xk ),  it  follows  that 

Sfc  =  E  [(afc_ i  o  (u-  xk-i)  +  Pyk )  (bk~ i  o(u-  xfc_i)  +  P^)T] 

=  E  [(afc-i  o  (u  -  xk-i))  (bk- i  o  (u  -  xk-i))T ]  +  E  [(afc_i  o  (u  -  xk-i))  ( Pzk)T ] 

+  E  [{Pyk)  {bk- i  o  (m  -  xfc_i))T]  +  E  [(Pt/fc)  (-P^fc)T] 

=  E  [(afc-i  o  (u  -  xk-i))  {bk-i  o  (u  —  xk-i))T]  +  E  [(afe_i  o  (u  -  Xfc_i))  {bk-i  o  a:fc_i)T]FrPr 
+  PFE  [(afc-i  o  Xfc_i)  {bk- 1  o  (M  -  Xfc_i))T]  +  PFE  [(afe_i  o  xk-i)  {bk- 1  o  aifc_i)T] FTPT 
=  (1  -  g)2  E  [afc-ift^J  +  g(l  -  g)  E  [diagK-i^)] 

+  g(l  -q)  E  [afc_i6fc_i]FTPT  -  g(l  -  g)  E  [diag(afe_iCi)]FrFT 
+  g(  1  -  g)PPE  [afc-ift^J  -  g(l  -  g)PFE  [diag(afc_i6^’_1 )] 

+  q2PF  E  [afc-i^_i]  F*PT  +  g(l  -  ?)EF  E  [diag(afc_i^_i)]PTPT.  (74) 

By  re-arranging  terms  in  (74)  and  observing  that  Sfc_i  =  E  [ak-ib^_)\  and  diag(Sfc_i)  =  E  [diag(afc_i6^_1 
the  result  in  (58)  follows.  Finally, 


Tfc+i  =  E  [yk+izk+1\  =  FE  [{ak  °  xk)(xk  o  bk)T]F 1 


=  F 


q-  E  [akbl]  +  g(l  -  g)  E  [diag(afc6j[)] 


=  F[g2Sfc  +  g(l  -  g)diag(Sfc)]Fr, 


(75) 


which  completes  the  proof.  ■ 

Although  omitted  in  the  statement  of  Lemma  3,  it  is  easy  to  see  that  the  dynamics  of  Ak  =  E [bkak] 
and  ©fc  =  E[zkyl]  can  also  be  obtained  by  noting  that  Ak  =  and  Qk  =  Tk. 
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V.  Convergence  Analysis  of  Robustified  Double-Iteration  Algorithm 

The  previous  Section  established  that  the  iterations  governing  the  evolution  of  T/-,  T /,.  and  5/,.  are 
identical  except  for  the  initial  conditions.  We  will  show  next  that  the  steady-state  solutions  of  these 
iterations  are  also  identical  up  to  a  multiplicative  constant.  To  see  this,  we  will  rewrite  (54),  (56), 
and  (58)  in  vector  form  using  Kronecker  products.  For  given  matrices  C,  A,  and  B  of  appropriate 
dimensions,  the  matrix  equation  C  =  AXB  (where  X  is  an  unknown  matrix)  can  be  rewritten  as  a  set 
of  linear  equations  of  the  form  (B  r  ®  A)x  =  c,  where  x  and  c  are  the  vectors  that  result  from  stacking 
the  columns  of  matrices  X  and  C  respectively,  and  ®  denotes  the  Kronecker  product4  of  matrices  [26]. 
Let  7k  be  the  vector  that  results  from  stacking  the  columns  of  Tk  and  %  the  vector  that  results  from 
stacking  the  columns  of  diag(Tfc).  Then,  it  can  be  easily  seen  that  (54)  can  be  rewritten  as 

7/,  =  [[qPF  +  (1  -  q)I]  ®  [qPF  +  (1  -  q)I]]  7*_i  +  [g(l  -<?)[/-  PF\  ®  [I  -  PF]] %_1}  k  >  1. 

(76) 

Let  G  be  a  diagonal  matrix  with  entries  G((Z  —  1  )n2  +  1,(1  —  1  )n2  + 1)  =  1,  V/  =  1,  2, ... ,  n2,  and  zero 
otherwise.  Then,  the  second  term  on  the  right  hand  side  of  (76)  can  be  written  as 

q(  1  -  q){[I  -  PF]  ®  [I  -  PF })%-!  =  q(l  -  q){[I  -  PF]  ®  [/  -  PF])Glk^,  k>  1,  (77) 

which  leads  us  to 

Ik  =  [[qPF  +  (1  -  q)I]  ®  [qPF  +  (1  -  q)I]  +  q{  1  -  q)([I  -  PF]  ®  [I  -  PF])G]lk_u  k  >  1.  (78) 

Let  ipk  and  (),  and  8k  be  the  vectors  that  result  from  stacking  the  columns  of  '1-7,  Ek  and  Ak  respectively. 
Then,  it  is  easy  to  see  that  the  same  recurrence  relation  as  in  (78)  governs  the  evolution  of  'ti+  and  £/,. . 

Theorem  1:  Let  P  G  Mnxn  be  a  column  stochastic  and  primitive  weight  matrix  associated  with  a 
directed  graph  Q  =  (V, 8},  with  V  =  {1,2, ...  ,n}  and  £  C  V  x  V.  Let  F  =  [In  In  ...  In]  G  Mnxn2, 
where  In  is  the  n  x  n  identity  matrix,  and  P  =  [EiPT  E2PT  . . .  EnPJ]T  G  Mn2xn,  where  each 
Ei  G  Rnxn,  i  G  {1,  2, . . .  ,n},  satisfies  Ei(i,i)  =  1  and  has  all  other  entries  equal  to  zero.  Then,  for 
any  q,  0  <  q  <  1,  the  matrix  II  defined  as 

fl  =  [qPF  +  (1  -  q)I\  ®  [qPF  +  (1  -  q)I\  +  q(  1  -  q)([I  -  PF]  ®  [/  -  PF])G  (79) 

is  column  stochastic,  and  it  has  a  single  eigenvalue  of  maximum  magnitude  at  value  one. 

Proof:  We  show  first  column  stochasticity  of  matrix  II.  Let  C  =  qPF  +  (1  —  q)I  and  D  =  I  —  PF, 
so  that  II  =  C +  q(l  —  q) (D® D)G.  We  will  establish  that  C®C  is  column  stochastic  and  also  show 

4The  Kronecker  product  of  matrices  A  —  [a,,]  £  Rmxn  an(j  B  =  [6,;,]  £  Rpx,J  is  defined  (see,  e.g.,  [26])  as  the  block  matrix 

ttii-B  .  . .  d\ nB 


dmiB  ...  a  mn  B 


A®  B  \= 


Xnq 


COORDINATED  SCIENCES  LABORATORY  TECHNICAL  REPORT  UILU-ENG-1 1-2207  (CRHC- 11-05) 


21 


that  the  column  sums  of  D®D  are  all  zero.  By  construction,  the  entries  of  the  ith  column  of  P  e  M"2xn 

are  all  zero,  with  the  possible  exception  of  the  ones  indexed  by  ((i  —  l)n  +  j,  i) ,  i,j  =  1,2 , ,n, 

2  ~ 

each  of  which  corresponds  to  the  (j,  i )  entry  of  matrix  P.  Then,  it  follows  that  Ya= i  =  Ej= 1  Pji  = 

1 ,Wi  =  1,2, ...  ,n2.  The  matrix  PF  e  M”2*™2  is  also  column  stochastic  by  construction,  as  it  results 

from  horizontally  concatenating  n  times  the  matrix  P,  i.e.,  PF  =  [P  P  ...  /  ];  therefore,  the  matrix  C 

is  also  column  stochastic.  The  kronecker  product  of  C  with  itself,  results  in  an  n4  x  n4  block  matrix  of 

the  form  C  <E>  C  =  [C\  C2  ■  ■  ■  Cn 2],  where  Cj  =  [c\jCT  c-^C'  . . .  cn2?C/  |  /  .  Then,  it  follows  that  the 

sum  of  the  entries  of  the  Ith  column  of  Cj  is  ^”=1  C3{m,  l)  =  (^''=1  c,:?)(^"=|  cr/).  Since  Y^=i  cv 

_  2  ~ 

and  £”=  M  cri  are  the  sum  of  the  entries  of  the  jth  and  Ith  columns  of  C  —  qPF  +  (1  —  q)Ini  (which  is 

4 

column  stochastic),  it  follows  that  El=i^i(m)0  =  ^  therefore,  C  is  also  column  stochastic. 

Since  PF  is  column  stochastic,  the  column-sums  of  D  =  /  —  PF  are  zero.  The  kronecker  product 
of  D  with  itself  is  of  the  form  D  ®  D  =  [Di  D2  . . .  Dn2],  where  Dj  =  [dijDT  d2jDT  . . .  dn2jDT]T. 

4  _  2  _  2 

Using  similar  arguments  as  above,  it  follows  that  Ylm=i  Dj  (' m,l )  =  (Er=idij)(Er=1dw)  =  0,  which 
implies  that  the  column- sums  of  D  0  D  are  zero.  The  only  thing  left  to  establish  that  II  is  column 
stochastic  is  to  show  that  all  entries  of  II  are  nonnegative  (from  where  it  immediately  follows  that 
II  =  C  0  C  +  q(  1  —  q){D  ®  D)G  is  column  stochastic).  We  argue  nonnegativity  of  II  as  follows: 
due  to  the  sparsity  structure  of  G  in  (79),  the  only  nonzero  entries  of  (D  <g)  D)G  will  be  in  columns 
(k  —  l)n2  +  k,  k  =  1,2, ... ,  n2;  thus  except  for  entries  in  these  columns,  the  entries  of  II  will  be  identical 
to  the  corresponding  entries  in  C ® C.  From  the  structure  of  PF,  entries  of  C and  q(l  —  q)(D ® D) 
can,  respectively,  take  one  of  the  following  three  forms: 

(qPij  +  (!  -  7))  (qpim  +  (1  -  ?)), 

9(1  -  q)(l-Pij)(l-Plm), 

or 

qpij  {qp  Im  +  (1-9)), 

-  9(1  -  q)PijO-~Plm), 

or 

q2PijPlm,  (84) 

9(1  -  q)pijPimi  (85) 

where  pt]  >  0  and  pirn  >  0  are  the  (i,j)  and  (l,m)  entries  of  matrix  P.  For  (80)  and  (81),  the 
corresponding  entry  of  IT  is  of  the  form 


(82) 

(83) 


(80) 

(81) 


(qpij  +  (1  -  9))  {qpim  +  (1  -  9))  +  9(1  -  9)(1  -  Pij){1  ~  Plm )  =  qPijPim  +  (1  -  9),  (86) 


COORDINATED  SCIENCES  LABORATORY  TECHNICAL  REPORT  UILU-ENG-1 1-2207  (CRHC- 11-05) 


22 


and  satisfies  0  <  qpijVim  I-  (1  —  q)  <  1  For  (82)  and  (83),  the  corresponding  entry  of  II  is  of  the  form 

qpij{qpim  +  (1  -  q))  -  q{  1  -  q)pij{  1  -  Pim)  =  qPijPim,  (87) 

and  satisfies  0  <  qpijPim  <  1-  For  (84)  and  (85),  the  corresponding  entry  of  II  is  of  the  form 

q2PijPim  +  q{  1  -  q)PijPim  =  qpijPim,  (88) 


and  satisfies  0  <  qp^pim  <  1. 

To  prove  the  second  assertion,  we  will  show  first  that  matrix  PF  can  be  written  via  a  permutation 
of  its  indices  in  the  form 


U  V 

0  w 


(89) 


where  U  is  an  irreducible  column  stochastic  matrix  and  lim^oo  Wk  =  0.  Since  PF  is  column  stochastic, 
we  can  assume  that  it  corresponds  to  the  weight  matrix  of  some  graph  Q  =  {V,£}.  We  will  show  that 
this  graph  has  a  single  recurrent  class  plus  a  few  transient  states,  from  which  the  decomposition  of  PF 
in  (89)  follows.  Let 


V  =  {(1, 1),  (2, 1), ... ,  (n,  1),  (1,  2),  (2,  2), . . . ,  (n,  2), . . . ,  (n,  n  —  1),  (1,  n),  (2,  n), . . . ,  (n,  n)}.  (90) 

From  the  structure  of  PF,  it  follows  that  for  any  node  (i,j)  E  V,  one-step  transitions  out  of  (i,j)  are 
to  nodes  of  the  form  (m,i),  with  i  E  M~,  where  A f~  is  the  set  in-neighbors  of  node  m  in  the  graph 
Q  (with  weight  matrix  P).  From  the  structure  of  PF,  it  also  follows  that  there  are  possibly  several 
rows  of  PF  with  all  entries  equal  to  zero,  which  means  that  a  node  (i,j)  that  is  associated  with  such 
row  cannot  be  reached  from  any  other  node;  however,  as  already  argued,  from  nodes  of  the  form 
it  is  possible  to  reach  nodes  of  the  form  (m,  i),  where  i  E  M~.  Clearly,  the  nodes  corresponding  to 
rows  with  all  entries  being  zero  are  transient.  Note  that  the  possibility  of  individual  nodes  that  cannot 
be  reached  from  any  other  node  being  disconnected  is  ruled  out  as  it  is  easy  to  see  the  only  nonzero 
diagonal  entries  of  PF  correspond  to  diagonal  entries  of  P,  which  are  strictly  smaller  than  one. 

Next  we  will  show  that  from  a  node  (i,j)  whose  corresponding  row  in  PF  has  some  nonzero  entries 
one  can  reach  any  other  node  (m,  l )  whose  corresponding  row  in  PF  has  some  nonzero  entries.  This 
means  that  all  non-transient  nodes  form  a  single  recurrent  class  (as  already  argued  all  nonzero  diagonal 
entries  are  strictly  smaller  than  one  which  means  there  cannot  be  absorbing  nodes).  This  follows  from 
the  fact  that  the  graph  Q  is  strongly  connected,  which  means  that  for  any  j,  l  E  V,  there  exists  a  path 
between  j  and  l.  Let  denote  the  nodes  traversed  along  the  path  between  j  and  /.  We  will 

show  next  that  for  any  two  non-transient  nodes  (i,j),  (r,  l)  E  V  there  exists  a  path.  As  already  argued, 
from  (i,j)  one  can  reach  in  a  single  hop  any  node  of  the  form  (m,  i),  where  m  is  a  neighbor  of  node  i 
in  the  graph  Q.  Since  i  \  is  the  first  node  traversed  in  the  path  between  j  and  l,  it  follows  that  (?'| .  i)  E  V 
can  be  reached  in  one  step  from  By  repeatedly  using  this  argument,  it  follows  that  the  sequence 
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of  nodes  {ii1  i),  ( i2 ,  ii), . . . ,  ( it ,  (r,  it)  forms  a  path  between  (i,j)  and  (r,  /),  which  means  that  any 

non-transient  node  can  be  reached  by  any  other  non-transient  node;  thus,  the  set  of  non-transient  nodes 
forms  a  single  recurrent  class.  Clearly,  the  vertex  set  V  can  be  decomposed  into  a  single  recurrent  class 
and  possibly  several  transient  nodes.  By  re-ordering  the  nodes,  it  follows  that  PF  can  be  rewritten  as 
in  (89)  (see,  e.g.,  [27,  p.  126]).  Furthermore,  since  Q  in  (89)  is  irreducible,  it  follows  that  qQ  +  (1  —  q)I 
(where  /  is  the  identity  matrix)  is  primitive.  It  follows  that  C  =  qPF  +  (1  —  q)  I  has  a  unique  largest 
eigenvalue  of  value  one,  i.e.,  Ai  =  1,  and  1  >  |A2|  >  •  ■  •  >  |An2|.  Let  a(C)  =  {Ai,  A2, . . . ,  An2}.  Then, 
a (C  0  C)  —  {AjAj,  i  —  1, . . .  ,n,  j  =  1, . . .  ,  n},  including  algebraic  multiplicities  in  both  cases  [26, 
p.  245].  Since  Ai  =  1  is  unique  (multiplicity  one)  and  | A*|  <1,  i  =  2, . . . ,  n2,  it  follows  that  the 
eigenvalue  of  C  0  C  =  [qPF  +  (1  —  q)I)  0  [ qPF  +  (1  —  q)I]  of  largest  magnitude  also  takes  value  1 
and  is  unique.  Since  C  0  C  is  column  stochastic  and  Ai  =  1  is  unique,  we  know  that  either  C  0  C  is 
also  primitive  or  it  can  be  decomposed  following  a  permutation  of  indices  to  the  form  [27,  p.  126]: 


L  M 
0  N 


(91) 


where  L  is  a  primitive  matrix  and  liui/.^oo  Nk  =  0. 

We  will  show  next  that  II  =  C  0  C  +  q(  1  —  q)(D  0  D)G  has  exactly  the  same  nonzero  entries  as 
C ®C  and  therefore  can  be  decomposed  following  the  same  permutation  of  indices  to  the  form  in  (91). 
As  argued  before,  due  to  the  sparsity  structure  of  G  in  (79),  the  only  nonzero  entries  of  (D  0  D)G  will 
be  in  columns  (k  —  1  )n2  +  k,  k  =  1,2,...,  n2,  thus  except  for  entries  in  the  aforementioned  columns, 
the  nonzero  entries  of  II  will  be  the  same  as  those  in  C  0  C.  For  all  other  columns  in  II  (that  include 
nonzero  entries  in  (D  0  D)G),  it  was  shown  in  (86)— (88)  that  the  nonzero  entries  of  II  are  strictly 
positive,  from  where  it  follows  that  II  has  the  same  sparsity  structure  as  C  0  C,  which  means  that  II 
can  also  be  decomposed  in  the  form  of  (91)  (for  some  matrices  L',  M',  N'),  and  the  resulting  upper-right 
block  is  also  a  primitive  matrix.  Therefore,  II  has  a  unique  largest  eigenvalue  at  one.  ■ 

The  following  two  lemmas  establish  that  the  first  and  second  moments  of  ap.  and  bk,  and  yk  and  zk 
converge  to  the  same  solution  up  to  a  scalar  multiplication.  These  two  lemmas  will  be  used  to  show 

e;=i  zoU) 


that  as  k  — >•  oo,  the  random  vector  vk  —  zk  —  oyk,  for  a  = 


E;=i  yo(j) 


— ,  will  converge  almost  surely  to 


v  =  0.  This  suggests  that,  as  k  ^  oo,  and  whenever  zk  is  nonzero,  each  node  i  can  obtain  an  estimate 
of  a  =  ^  lZyv(j)  by  calculating  the  ratio  yk(i) / zk(i).  We  will  also  show  that,  in  fact,  zk  will  be  larger 
than  some  threshold  infinitely  often. 


Lemma  4:  The  first  moments  of  ak  and  bk  (also  yk  and  zk  asymptotically  converge  to  the  same 
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solution  up  to  scalar  multiplication: 


lim  zk  =  ot  lim  yk, 

h—too  k — yoo 

lim  bk  =  a  lim  ak, 

k — ^oo  k — yoo 


(92) 

(93) 


where  a  = 


EJ=i*o(j) 

EUvour 


Proof:  In  Lemma  1,  it  was  shown  that  yk+1  =  [qP  +  (l  —  q)l]yk  and  zk+i  =  [qP  +  (l—q)l]zk  with 
Vi  —  QVo,  and  z\  =  qz0.  Since  P  is  column  stochastic  and  primitive,  it  follows  that  [qP+(l  —  q)l]  is  also 
column  stochastic  and  primitive.  Thus,  lim^oo  zk  =  a  lini/.^00  where  from  the  column  stochasticity 
property  it  follows  that  =  ?(E"=i  zo(j))  and  E”=i  Vkij)  =  q(  E”=i  Vq(J))  >  >  !;  this 

zn(Y) 

implies  that  a  =  ~  ,  which  establishes  (92). 

r  Ei=iyou)  _  _ 

By  noting  that  a0  =  Pyo,  and  b0  =  Py0,  and  using  the  fact  that  P  is  column  stochastic,  it  follows 

that  E;=i«*(i)  =  e;=i«o(j)  =  E;=i?/o(i)  and  e;=i^(j)  =  e;=Jo(j)  =  euz^)-  since 

gPF  +  (1  —  q)I  (i.e.,  the  matrix  that  governs  the  dynamics  of  ak  and  bk)  is  column  stochastic  and,  as 
shown  in  the  proof  of  Theorem  1,  has  a  single  largest  eigenvalue  at  value  1,  a  similar  development  to 
the  one  above  can  be  used  to  show  (93).  ■ 

Lemma  5:  Define  wk  =  bk  —  aak  and  denote  by  Xk  the  vector  that  results  from  stacking  the  columns 
of  Xk  :  =  E [wkwk],  Then,  it  follows  that 


Xk  =  Hxk-1,  (94) 

4 

with  xo  =  -0o  +  «27o  -  «(£o  +  4)  and  E"=i  Xo(Q  =  0. 

Proof:  Since  Xk  :=  E [wkwk]  =  E [bkbl]  +  a2E[akal]  -  a(E[akbk]  +  E [bkak])  =  \Dfe  +  a2Tk  - 
a(Ek  +  Afc),  it  follows  that  Xk  =  f>k  +  oi2^k  —  a(fk  +  5k).  From  (78)  and  subsequent  discussion,  it 
follows  that  7fc  =  n7fc_i,  0fc  =  II0fc_i,  =  nffc_i,  and  6k  =  IL5fc-i,  thus  Xk  =  Rfk-i  +  a2n7fc_i  ~ 

a(nefc-i  +  n<$fc-i)  =  n(0fc_i  +  «27fc_i  -  a(ffc_  i  +  4- 1))  =  n^fc-i. 
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In  Lemma  3,  it  was  shown  that  r0  =  Py0yQpT,  'ho  =  PzqZqPt ,  and  H0  =  Py0ZQpT  =  Aq  .  Since 
7o>  Po,  £o>  and  5o  result  from  stacking  the  columns  of  r0,  'ho,  50,  and  A0,  it  follows  that 

ro(*,j)  =  (5^Z/o(®) 

1=1  i=  1  j= 1  \  i=l 

^0  (i,j)  =  W^z0(i 

1=1  2=1  j= 1  \  2=1 

n4  n2  n?  /  n  \  /  n 

=  ^2^2^o(hj)=  (5Zj/°(*)] 

1=1  i=l  j=l  \  i=l  )  \i=l 

n4  n?  n2  /  n  \  /  n  \ 

)  (X^o(X>)’  (97) 

1=1  i= 1  j=l  \i=l  /  \i= 1  / 

where  the  last  equality  is  obtained  by  taking  into  account  that  i)  matrix  P  is  column  stochastic  by 
construction,  and  ii)  for  any  a,  b  e  Mn,  we  have  that  E”=i  E™= i  abT(i,j)  —  (YPi= i  a*)(ELi  Since 

a  =  follows  that  E"=i  AA)(0  =  EI=i(^o(0  +  «27o(0  -  a(fo(0  +  <U0))  =  °-  ■ 

Theorem  2:  Let  yk  and  Zk  be  the  random  vectors  that  result  from  iterations  (28)-(29)  and  (30)— (3 1). 
Define  vk  =  Zk  —  ayk,  where  a  =  (77 1  yyP  Then,  ||iy  ||oc  — >  0  almost  surely.  Furthermore,  for  every 
j,  vk{Jj)  — y  0  as  k  — y  oo  almost  surely  (i.e.,  for  every  j,  lim/;_>.00  vk(j)  =  0  with  probability  one). 

Proof:  The  result  follows  from  the  first  Borel-Cantelli  lemma  [28,  Theorem  7.3.10].  For  all  k  >  0 
and  all  e  >  0,  define  the  event  Ek{e)  =  {||ffc||oo  >  e}.  We  will  first  establish  an  upper  bound  on 
EtioPH-M6)}  by  noting  that  Pr{£t(e)}  =  Pr{ ||t>t||oo  >  «}  <  E'1";1  b  thus  E7o Pr{-Ei(e)}  < 
1E£oE[IKIU  <  fEr=oE[IKIIJ.  Note  that  E [|k||2]  =  (E[u>t]))‘/2  =  (trace(E[urf]))V2  = 

(trace(E[^fc^J])  +  a2trace(E[i/fc|/J])  —  2atrace(E[yfc^])1'/2.  We  will  next  show  that  E[||ufc||2]  — *  0  as 
k  — *  oo  geometrically  fast.  Using  Lemma  3,  it  can  be  established  that  E [vkvk\  =  E [zkzk]+a2  E [ykyk\  — 
a{E[ykz% }  +  E[zkyl]))  =  +  q(  1  -  g)diag(Xfc_i)]Fr  where  Xk_x  =  E [wkwf]  as  defined  in 

Lemma  5,  thus  the  evolution  of  E[vkvk]  is  governed  by  the  evolution  of  Xk_i  or  by  Xk-i  (the  vector  that 
results  from  stacking  the  columns  of  Xk_f).  In  Theorem  1,  we  showed  that  If  has  a  unique  eigenvector 
(with  all  entries  strictly  positive)  associated  to  the  largest  eigenvalue  Ai  =  1.  Then,  the  solution  of  (94)  is 
unique  and  equal  to  this  eigenvector  (up  to  scalar  multiplication).  Since  II  is  a  column  stochastic  matrix, 

4  _  4 

and  Lemma  5  established  that  E"=  iXo(0  =  0,  it  follows  that  ff'l-A  XkO)  =  0.  k  >  0,  and  therefore 
lim k^ooXkP)  =  0,  V7.  Additionally,  it  is  well-known  that  the  convergence  of  (94)  is  geometric  with  a 
rate  of  convergence  given  by  the  eigenvalue  A2  of  II  with  the  second  largest  modulus,  which  satisfies 
| A2 1  <  Ai  =  1  (see,  e.g.,  [27]).  Thus,  we  have  established  that  Xk( 0  —3 ►  0,  VZ,  geometrically  fast, 
from  where  it  follows  that  all  the  entries  of  E[vkvk]  go  to  zero  also  geometrically  fast.  Therefore,  the 
trace (E[vkvk])  also  goes  to  zero  geometrically  fast,  so  that  E  [||vfc||2]  also  goes  to  0  geometrically  fast. 
It  immediately  follows  that  Efclo  ®  [Ikfclb]  <  00  and  therefore  E/T=o  Pf{  Halloo  >  e}  <  oo.  Then, 
from  the  first  Borel-Cantelli  lemma  PrdltifcUoo  >  e  infinitely  often}  =  0  (or  Pr{||'Ufc||00  >  e  i.o.}  =  0). 


Z)7o(0  = 
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Finally,  since,  for  every  j.  Halloo  >  \vk(j)\,  then,  for  every  j,  PrdlufcHoo  >  e}  >  Pr{|nfc(j)|  >  e},  and 
thus,  for  every  j,  F>r{|t’fc(i)l  >  e}  <  oo.  Then,  by  Theorem  7.2.4.C  of  [28],  for  every  j,  vk(j)  — >  0 
almost  surely.  ■ 

Theorem  2  has  established  that,  in  the  limit  as  the  number  of  iterations  k  becomes  large,  the  values 
of  vectors  yk  and  zk  will  be  perfectly  aligned  so  that  zk  —  ayk  =  0  with  probability  one.  Thus,  in  this 
limiting  case,  each  node  j  can  calculate  the  value  of  P  by  taking  the  ratio  as  long  as  zk(j )  f  0. 
Note  that,  as  also  evidenced  by  the  simulations  provided  for  the  small  network  of  Fig.  1  (e.g.,  the  plots 
on  the  left  and  in  the  middle  for  Figure  3),  the  vectors  yk  and  zk  do  not  converge  in  any  way;5  however, 
the  values  yk  and  zk  become  perfectly  aligned  (with  probability  one),  allowing  each  node  j  to  calculate 
«  =  TTjj-  The  only  problem  here  arises  when  yk(j)  and  zk  {] )  have  both  value  zero,  which  does  not 
constitute  a  violation  of  zk  —  ayk  =  0,  but  clearly  does  not  allow  node  j  to  calculate  the  desired  value 
4.  This  is  evidenced  also  in  the  simulations  provided  for  the  small  network  of  Fig.  1:  for  example,  in 
the  plots  in  Fig.  4,  the  values  of  yk(j )  and  zk(j)  often  go  to  zero  (simultaneously)  leaving  their  ratio 
undefined.6  The  next  two  theorems  essentially  establish  that  zk  (j ) ,  j  =  1,2,....  n,  will  be  greater  than 
zero  (in  fact,  greater  than  a  constant  C  that  will  be  specified)  infinitely  often.  Note  that,  in  subsequent 
developments,  zk(j )  is  denoted  with  z3 [k]  in  order  to  remain  close  to  the  notation  in  (15)— (17). 

Theorem  3:  Consider  a  (possibly  directed)  strongly  connected  graph  Q  =  (V,  £)  and  the  iteration 
in  (15)— (17),  where  x]t [/;■] ,  (j,  i)  £  8,  k  =  0,1,2,...,  are  independent  identically  distributed  (i.i.d.) 
indicator  R.V.’s  as  defined  in  (1),  i.e.,  x3l[k]  =  1  with  probability  q  and  x]t [k]  =  0  with  probability 
1  —  q,  independently  between  (),  i)  £  8  and  independently  for  different  k.  For  every  j  =  1,2, ... ,  n, 
define  the  event  EJk  =  {zj[kn\  >  C},  k  >  1,  where  C  =  {n+m)^+ax)n-i .  =  maxiGV{£>+}, 

n  =  | V | ,  and  m  =  \8\.  Let  Q3k  denote  the  indicator  of  the  event  Ej,.  k  >  1,  i.e.,  Q3k  =  1  whenever 
El  k>  1  occurs,  and  Q  =  0  otherwise.  Then,  whatever  Ci,  C2,  •  •  • ,  Cfe-i>  we  have  that 

Pr {Zj[{k  +  1  )n]  >  C  |  Cl  Ci_!,  ■■■,£}>  Qn ,  Vj.  (98) 

Proof:  Note  that  the  iteration  in  (15)  to  (17)  involves  nonnegative  quantities:  since  for  every  j, 
Zj[ 0]  >  0,  Vj,  it  follows  from  (30)-(31)  that,  for  every  j,  Zj[k\  >0,  k  >  0.  Then,  it  is  not  hard  to 
establish  that  the  total  mass7  M.k+i  in  the  system,  defined  as 

n 

Mk+ 1  :=  yjZj[k  +  1]  +  y  (Cji[k]  -  Tji[k  -  1])(1  -  Xji[k])  ,  (99) 

i=1 

5Earlier,  we  established  that,  for  large  k.  the  quantities  E[yk],  E[zk],  E[ykUk\,  E[zkZk]  and  E[ykZk]  converge,  but  this  does  not  imply 
any  convergence  for  the  values  of  yk  or  zk- 

6Since  in  the  simulations  for  the  plots  in  Fig.  4,  each  packet  (including  self-packets)  can  be  dropped  with  probability  1  —  q  at  iteration 
k ,  there  is  a  nonzero  probability  that  all  packets  destined  for  node  j  will  be  dropped,  causing  both  of  its  values  at  the  next  iteration 
(yk+i{j)  and  Zk+i(j))  to  be  zero.  For  instance,  in  the  simulation  of  Fig.  4,  2^(1)  will  be  zero  with  probability  at  least  (1  —  q)2  =  0.81 
because  node  1  will  have  value  zero  if  both  packets  destined  for  it  (including  the  self-packet)  are  dropped. 

7This  notion  is  discussed  in  great  detail  in  Part  II  of  this  paper. 
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satisfies 


Mk+i  —  n  ,  for  all  k  —  0, 1,  2, ...  . 


[This  follows  from  the  fact  that  M0  =  V "  ,  2y[0]  =  n  and  the  observation  that 


Mk+i  ■=  X  zi  +  X  fo*  lkl  -  Ta  lk  ~  !] )  ( : 1  -  xa  [ k ] ) 

i=1 

=  X  (vjM  -  rAk  -  l\)xn[k]  +  X  (aAk\ -  -  1i)(1  - x. 

=  Mk)  -  rAk  - 


j* 


X  (  a3i  [k  ~ 


+  -  aAk  -  l]x3i[k  -  1]  -  Tji[k  -  2](1  -  Xji[k  -  1]) 


X^'M  +  X  Aji\k  -  !]  -  Tjilk  -  2])(i  -  Xji[k  -  1])  , 
j=1  0,0  g£ 


which  is  equal  to  M.k .] 

The  definition  of  M.k+i  in  (99)  involves  the  summation  of  n  +  rri  nonnegative  quantities,  namely, 
Zj[k  +  1]  for  j  =  1,  2, n  and  m,ji[k  +  1]  :=  (crjjfc]  —  r^/c  —  1] ) (1  —  XjJ/c])  for  (j,i)  e  We  can  think 
of  these  quantities  as  follows:  Zj[k  +  1]  is  the  mass  at  node  j,  whereas  rriji[k  +  1]  is  the  mass  waiting  to 
get  transferred  to  node  j  from  node  i.  Since  all  of  these  quantities  are  nonnegative,  at  least  one  of  them 
is  larger  or  equal  to  Regardless  of  whether  this  quantity  is  associated  with  a  node  (say  node  j*) 
or  a  link  (say  link  (j*,  i*)),  this  mass  has  at  least  one  way  of  reaching  any  node  i  of  interest  in  graph 
Q  via  a  path  of  length  at  most  n  —  1  (because  the  graph  Q  is  strongly  connected):  in  particular,  there 
is  at  least  one  path  of  length  at  most  n  —  1  from  node  j*  to  node  i  and  all  the  links  in  this  path  have 
weight  at  least  —2 — .  If  all  these  links  are  activated,  which  occurs  with  probability  qn~l  (qn  in  the  case 

^max 

of  link  ( j*,i *)  because  the  mass  needs  to  first  transfer  to  j*),  then  a  fraction  ( — ^ — ) r'~ 1  of  the  mass 

l^max 

will  transfer  to  node  i  in  at  most  n  steps.  Then,  since  for  every  j,  Zj[k]  >0,  k  >  0,  independently  of 
the  values  of  Zj[ln\,  l  =  1,  2, . . . ,  k,  Pr {zj[(k  +  1  )n\  >  C  \  (Jk,  Ci-n  ■  ■  •  >  Ci }  >  (f'  obtains,  whatever 
Ci,  C2,  -  -  - ,  Cfc-1-  Finally,  for  every  j  ,  Pr {zj[(k  +  l)n]  >  C  \  Cl  (Lv  ■  ■  ■  >  (1}  =  1  “  p *{zj[(k  +  1)n]  < 
C  |  (I,  Cl-i,  ■■■,({}<  l  —  Pr{zj[(k+l)n\  =  0  |  (Jk,  (Jk_v  ■  ■ . ,  Cl}  <  1  ~qVj  ,  where  Vj  is  the  in-degree 
of  node  j.  ■ 

Given  a  sequence  of  events  E\ ,  E2, . . . ,  En, . . .  defined  on  some  probability  space,  the  next  theorem 
(which  we  do  not  prove)  states  the  1912  Borel  criterion  for  establishing  whether  the  event  that  infinitely 
many  of  the  Ek  occur,  denoted  by  {Ek  i.o},  will  occur  with  probability  one  or  zero  (see,  e.g.,  [29], 
[30]).  This  result,  together  with  the  result  in  Theorem  3  will  be  used  to  establish  that,  for  every  j,  the 
event  E:k  =  {zj[kn\  >  C},  k  >  1  occurs  infinitely  often. 

Theorem  4:  Let  {Ek},  k  =  1,  2, . . . ,  be  a  sequence  of  events  defined  on  some  probability  space.  Let 
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Qk  be  the  indicator  function  of  the  event  Ek.  Let  Pr{C/,+)  |  Qk,  ,  Ci}  denote  the  conditional  prob¬ 

ability  of  the  event  Ek+l  given  the  outcome  of  previous  trials.  If  0  <  V'k  <  Pr{^fc+i  I  C k,  Cfc-i,  ■  ■  ■ ,  Ci}  < 
p"  for  every  k,  whatever  Ci,  C2,  ■  ■  ■ ,  (k,  then  i)  Pr{Efc  i.o.}  =  0  if  Y^kLiPk  —  00 >  and  h)  Pr {Ek  i.o.}  =  1 

tfE£iK  =  °°. 

Theorem  5:  Consider  a  (possibly  directed)  strongly  connected  graph  Q  =  (  V.  S)  and  the  iteration 
in  (15)— (17).  For  every  j  =  1,2,  define  the  event  EJk  =  {zj[kn\  >  C},  k  >  1,  where  C  = 

(n+m)(p+ax)"-i >  pmax  =  maxieV{X>t},  n  =  |V|,  and  m  =  \£\.  Then,  Pr {Ek  i.o.}  =  1. 

Proof:  Theorem  3  established  that,  for  every  j,  Pr {zj[(k  +  l)n]  >  C  \  Cl,  Cl-n  ■  ■  ■  ■>  C{}  —  'f  - 
Define  p'k  =  qn,  then  it  follows  that  YlkLiPk  =  00 >  and  by  the  second  assertion  in  Theorem  4,  we 
conclude  that,  for  every  j,  Pr {E3k  i.o.}  =  1.  ■ 

The  final  piece  is  to  establish  that  whenever  Zj  [k]  >  C,  which  occurs  infinitely  often,  each  node  will 
be  able  to  calculate  an  estimate  off;  by  calculating  the  ratio  ]jj[k]/zj[k]  and  this  estimate  will  converge 
to  1/a  as  k  goes  to  infinity. 

Theorem  6:  For  each  j,  let  k  —  ti,  f2,  ■  ■  •  be  an  increase  sequence  of  time  steps  for  which  Zj[k]  >  C. 
Then,  almost  surely 

Vj\Pn\_  _ 

J  a 


lim 

n— >00 


Zj  [tr 


=  0. 


(100) 


Proof:  Since  Zj  [k]  >  C  for  k  —  t\,  t2, . . . ,  it  follows  that 


Uj  [t"n 


Zj  [tn  ] 


< 


^ .  Also,  in  the  proof 


of  Theorem  3,  we  established  that  A4k  —  n,  k  >  0,  from  where  it  follows  that  Zj[tn\  <  n,  therefore 
vzJ |*"j  —  4  >  _  jn  Theorem  2,  we  established  that  ay3 [k]  —  Zj[k]\  — >  0  almost  surely,  which 

implies  that  the  subsequence  \ayj[tn\  —  z3 [tn]  — y  0  almost  surely,  then  since  C  <  n,  we  have  that 


lim 

Vj  [tn]  1 

<  lim 

UVj  \pn\  Zj  [tn\ 

n— >oo 

Zj[tn\  a 

71— »•  OO 

aC 

=  0 


(101) 


almost  surely. 


VI.  Concluding  Remarks 

In  this  paper,  we  proposed  a  method  to  ensure  robustness  of  a  class  of  linear-iterative  distributed 
algorithms  against  unreliable  communication  links  that  may  drop  packets.  We  used  statistical-moment 
analysis  and  the  Borel-Cantelli  lemmas  to  establish  the  correctness  of  the  proposed  robustified  algorithm. 
In  Part  II  of  this  paper,  we  establish  similar  convergence  properties  by  recasting  the  problem  as  a  finite 
inhomogeneous  Markov  and  using  coefficients  of  ergodicity  commonly  to  used  in  analyzing  this  type 
of  Markov  chains. 
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